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An effective scheme within two displaced bosonic operators with equal positive and negative 
displacements is extended to study qubit-oscillator systems analytically in an unified way. Many 
previous analytical treatments, such as generalized rotating-wave approximation (GRWA) [Phys. 
Rev. Lett. 99, 173601 (2007)] and an expansion in the qubit tunneling matrix element in the deep 
f*^ ' strong coupling regime [Phys. Rev. Lett. 105, 263603 (2010)] can be recovered straightforwardly 

i within the present scheme. Moreover, further improving GRWA and extension to the finite-bias case 

are implemented easily. The analytical expressions are then derived explicitly and uniquely, which 
work well in a wide range of the coupling strengthes, detunings, and static bias including the recent 
experimentally accessible parameters. 
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I. INTRODUCTION 



Matter-matter interaction is fundamental and ubiquitous in modern physics ranging from quantum optics, quantum 
information science to condensed matter physics. The simplest paradigm is a two-level atom (qubit) coupled to the 
electromagnetic mode of a cavity (oscillator). In the strong coupling regime where the coupling strength g/ui (lu is 
, ^ the cavity frequency) between the atom and the cavity mode exceeds the loss rates, the atom and the cavity can 
repeatedly exchange excitations before coherence is lost. The Rabi oscillations can be observed in this strong coupling 
£SJ | atom-cavity system, which is usually called as cavity quantum electrodynamics (QED) Typically, the coupling 
. strength in cavity QED reaches g/uj ~ 10~ 6 . It can be described by the well-known Jaynes-Cummings (JC) model Q 
C*~) ' without the rotating-wave approximation (RWA). 

Recently, for superconducting qubits, a one-dimensional (ID) transmission line resonator Q or a LC circuit 0-0| 
can play a role of the cavity, which is known today as circuit QED. More recently, LC resonator inductively coupled to 
a superconducting qubit [8l4lOl] has been realized experimentally. The qubit-resonator coupling has been strengthened 
from g/ui ~ 1CP 3 in the earlier realization Q, a few percentage later [1, 0] , to most recent ten percentages [8l-fl(|. 
Due to the ultra-strong coupling strength g/oj ~ 0.1, evidence for the breakdown of the RWA has been provid ed M . 
Recently, some works have been devoted to this qubit-oscillator system in the ultra-strong coupling regime [ 1 lM 1 5j j . 

Actually, the JC model without the RWA in a wide coupling regime has been studied extensively for more than 
40 years. By polaronic-like transformations or displaced operators, various analytical and numerical approaches have 
been developed in recent years, an incomplete list is given by Refs | 1 6j -|23ll - Very accurate or exact solutions have been 
obtained. Various forms of RWA energies have been developed [2J, [25| and extensions to the N-level case have been 
performed recently [26]. Most works have been mainly devoted to the zero static bias. The analytical expression with 
high accuracy for the qubit-Oscillator systems with both zero and finite static bias should be of practical interest. 

In this paper, by using two displaced bosonic operators with equal positive and negative displacements, we can 
recover many previous analytical treatments for both zero and finite static bias within the same scheme. What is 
more, we can extend analytically the previous generalized RWA (GRWA) to the finite static bias. Beyond the GRWA 
for the zero bias is also performed. The expression is uniquely given and the results are very close to exact ones for 
wide range of the model parameters which cover the present-day experimentally accessible parameters. 



II. MODEL AND EXACT SOLUTION 



The Hamiltonian for a superconducting qubit coupled to a harmonic oscillator in circuit QED consists of three 
parts [H [3. The first one is the interaction between the qubit and the LC resonator, which is described by 

Hint = %(a t + a)a z , (1) 
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where a\ a are the photon creation and annihilation operators in the basis of Fock states of the LC resonator, g is 
the qubit-cavity coupling constant. The RWA has not been employed here. The effective Hamiltonian for the qubit 
can be written as the standard one for a two-level system 

H = --(ea z + Aa x ), (2) 

where e and A are qubit static bias and tunneling matrix element. In the recent circuit QED || 0] operating in 
the ultra-strong coupling regime, they describe the transition frequency of the flux qubit and the tunneling coupling 
between the two persistent current states, e = I p ($ — <&o/2) with I p the persistent current in the qubit loop, $ an 
externally applied magnetic flux, and $o the flux quantum. In contrast to atomic cavity QED systems, e is easily 
tunable in circuit QED systems using superconducting qubit. In the above two equations, the Pauli matrix notations 
(Jk{k = x, y, z) are used in the basis of the two persistent current states. The third one is LC resonator uia^a with 
single mode frequency uj. Then the Hamiltonian for the whole system reads (H = lj = 1) 

H = -^(ea z + Aa x )+a^a + g(a? + a)a z . (3) 



Motivated by the work in the Dicke model 17[ , we have introduced two displaced bosonic operators with equal 
positive and negative displacements in this system fl5l| 

A = a + g, B = a - g, (4) 

then the Hamiltonian can be written in the following matrix form 

H-( A " A ~9 2 - £ 1 -A/2 \ () 
U -{ -A/2 Bt B _ ff 2 + 1 j ■ W 

Note that the linear term for the original bosonic operator a'(a) is removed, and only the number operators A + A 
and B + B are left. Therefore the wavefunction can be expanded in terms of these new operators as 



En=o(~ l ) n dn\n) B , 



(6) 



where N tr is the truanted number. For A operator, we have 

vn! vni 
\0)a = ^ 2 - sat) |0) a , (8) 
B operator has the same properties. Inserting Eqs. (6) and (7) into the Schro dinger equation, we have 

(m - g 2 - |j c m - D mn d n = Ec m , (9) 

n=0 

(m - g 2 + |) d m ~ ^ DmnCn = Ed m , (10) 
ra=0 

where 

D mn - j(-l) m B(m\\n) A , (11) 

B {m\\n) A = (2 3 )»- m exp(-2 ff 2 )^i— (V), (12) 

for n > m, L^ _m (a;) is Laguerre polynomial, D mn = D nm . 

Based on Eqs. (j9j and (|10j) . we have given numerically exact solutions to the qubit-oscillator system with any 
finite static bias e[15j. In this paper, we alternatively present some analytical results in the framework of the above 
formalism. One can see that some recent analytical results by other authors are explicitly covered in the present 
framework. Moreover, the present scheme is more convenient to perform further analytical studies. 
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III. ANALYTICAL TREATMENTS 
A. Variational study for e = 

To have a sense of two displaced bosonic operators Eq. (j4]) , we relax the displacement to be a variational parameter 

a, 

A = a + a,B = a — a, (13) 

then study the unbiased Hamiltonian (e = 0) variationally. Suppose that the trial state is the vacuum state in these 
displaced operators as the following 



l*o> = 



\0)a 



The energy expectation is derived as 



Minimizing the energy gives 



E = -ye- 2 " 2 -2ga + a 2 . (14) 



Aae- 2a -g + a = 0. (15) 
In the weak coupling limit, we can obtain the variational parameter and the ground state (GS) energy respectively 

a = — ~~7") (16) 
1 + A v ' 



A 

E = - — exp 



-2 



2 
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1 + A 



2 1 + 2A 

9 5, (17) 

(1 + A) 2 ' 



which are exactly the same as Eqs. (7) and (8) obtained in Ref. [27] . 

In the strong-coupling limit, the first term in Eq. (|14p. which is originated from the qubit tunneling, is too small 
and can be neglected, then we simply have 

a = g, (18) 

and the GS energy 

E SCL = -g 2 . (19) 



For the arbitrary coupling, one can solve Eq. (fTSj) consistently and the reasonable GS energy will be derived, which 
is not shown here. 

B. Perturbation theory based on the exact solution in the strong coupling limit 

Note above that in the strong coupling limit, the variational parameter is just exactly a = g. It can be also readily 
obtained by neglecting the qubit tunneling term —^Aa x in Hamiltonian ([3]) with zero static-bias. In this case, based 
on the displaced operators A and B, the eigenstates are easily obtained as 



\(o) _ ( \m) 

±(-ir\m) B 



(20) 



and the eigenvalues are — m — g 2 for the m state. Note that the eigenstates are twofold degenerate. 

Next, considering H' — —^Aa x as a perturbation, within the second-order perturbation theory, we can readily 
derive the eigenenergy with even (odd) parity for zero qubit static bias as 



E% = m- g 2 + (-l) m D mm + £ 1±^L. (21) 



m — n 
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FIG. 1: The energy levels for zero bias e = at (a) g = 0.1, (b) g = 0.5 , (c) g = 1.0, and (d) g = 1.5. The present BRWA 
results (black filled circles) are compared to the exact (solid lines), DSC(red dashed lines), and GRWA (blue dashed lines) ones. 



It is interesting to note that it is just the same as Eq. (5) in Ref. [21| obtained by Casanova et al. in the deep strong 
coupling (DSC) regime of the JC model. 

The energy levels by Eq. (|2ip against the qubit-oscillator detuning 6 = A — 1 for g = 0.1, 0.5, 1.0, and 1.5, ranging 
from weak to deep strong coupling, are displayed in Fig. [TJ The numerically exact results from Eq. (|9]) and (|10p 
are also collected as a benchmark. It is found that the DSC results are especially suited to the DSC regime or small 
detunings. Note that Casanova et al just focused on the investigation in the DSC regime {g/uJ — 2) or small detunings 
(A < 0.5). At weak coupling g — 0.1, it is shown in Fig. Q] (a) that, even for the negative detuning S, the DSC 
deviates from the exact ones. However, in the present experimentally accessible systems, the maximum value for the 
coupling strength is generally realized in the superconducting flux qubit coupled to a circuit resonant which is 
only around g = 0.1, to our knowledge. So it should be practically interesting to find a good solution in this coupling 
regime. 

For any value of the qubit bias e, the Hamiltonian ([5]) with a vanishing tunneling element A = can be diagonalizcd 
in terms of two eigenstates |f, m) A and ||, m) B with |f ) (||)) the eigenstate of a z , the corresponding eigenvalues are 



E 



(0) 
tA,i 



1 

±—e + m 



For finite A, the perturbative matrix elements becomes 

- 77 b (l,m\Aa x \f,n) A 



(-l) m D 



(22) 



(23) 



Note that these two euqations are exactly the same as Eqs. (7), (8) in Ref. |14 |. Then the full Hamiltonian can be 
diagonalized perturbatively to second-order in A by using Van Vleck perturbation (VVP) theory as outlined in Ref. 
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[Tij ]. The eigenvalue is given by 



I 1 

#m = m + - - g 2 + - 



n 2 



D 2 
u nk 



2 ^— ' V £ + f?7, — e + k — n 



1 



T 2 



\ 



/c=0,fc^m±Z 

m + l(l> 0),m = 0, 1,2,. 



r>2 



e + m— k e + k — n 



AD 2 



(24) 



which is the same as Eqs . (12) for VVP in Ref. [14j . So the VVP for finite bias can be also recovered easily in the 
present scheme. It can be reduced to the zero-bias case Eq. (|2"Tj) by set e = (m = n). 

It has been shown [l4| that VVP works very well in the deep strong coupling or large static bias. It is consistent 
with the fact that the unperturbative Hamiltonian includes the qubit-oscillator interaction and qubit bias. What 
happen for the accessible parameters of the present-day experiments? In addition, VVP at small static bias e < 1 has 
not been discussed either so far, which might however be more important. 

Here, we calculate energy levels in the VVP for different static bias s < 1, which are exhibited in Fig. [5] Compare to 
the exact ones, one can find that VVP deviates strongly with the increase of the tunneling parameters A in the a wide 
coupling regime g < 0.5, and become more pronounced at small static bias. Especially, around the experimentally 
accessible coupling strength around g — 0.1, VVP becomes worse considerably. Therefore a new analytical treatment 
is highly desirable. 



C. Analytical approximations at different levels 



In the framework of Eqs. ([9]) and (|10[) . analytical approximations can be performed systematically. First, as a 
zero-order approximation (ZOA), we omit the off-diagonal terms and have 



.9" 



5"* 



m 



E dr, 



0, 
0. 



Nonzero coefficients will give the following equation 

/ 9 e 



-E 



)( 



The eigenvalues are then given by 



The corresponding eigenstate is 



E+ = m - 



m — g 



9 T 2 



2~ E l D m,m - °' 



AD 2 

m.mi 



(-l) m D 
(m — g 2 - 



\m) 



#4 



(25) 



(26) 



The ZOA energies with zero static bias are just the three terms obtained in Eq. (|2"Tj) . In Fig. [21 we also plot the 
ZOA energy levels against the coupling constant g for several static bias. It is demonstrated from the upper and 
middle panel that for small static bias ( e ^ 0.5), ZOA is almost equivalent to the VVP in all parameters. If the high 
accuracy is not required, the simple expression of the eigensolutions in the ZOA should be practically very useful, at 
least as a preliminary estimate of some physical quantities . 

The approximation can be easily improved step by step with the consideration of more off-diagonal ele- 
ments in the present formalism. The first-order approximation (FOA) is performed by selecting two coefficients 
Cm d m c m +i d m+ i . The determinants for any m is given by 



n-(E) -D mm 

-An+l,m 0~ +1 (i£) 

— -Dm+l,m -^m+l,m+l 





— Dm+l.m+l 



= 0, 



(27) 



6 



where 



n- m (E) = m-g 2 -E-^ Sl+(E) = m-g 2 -E+ £ -. 



(28) 



Some roots of this quartic equation will give the energy levels. The analytical expression might be a little bit complicate 
but should be given unambiguously. 

We first revisit the zero-bias case e = 0. In this case, due to the parity symmetry, we can set d n = ±c„, then both 
equations give (m — g 2 ) c m =p X) 
form 



n =o D mn c n = Ec m . In the FOA, the determinant takes the following 2-by-2 block 



(m - g 2 - E =f An,™) TAn,m+i 
TD m +i, m (m + 1 - £ 



EtD 



m+l,m+l « 



= o, 



where the sign — (+) is for even (odd) parity. We can readily have two roots for even parity 

Eg>*> = m-(/ 2 + i- i (An+l.m+l + An.m) ± ^ \j [1 + (An.m - £> m+ l, m+ l)] 2 + 4£>^ m+1 , 



(29) 



and other two roots for odd parity 

E^ = m -g 2 + ± + ± (An+i.m+i + D m>m ) ± ^[1 - (£> m>m - £> m+ i,m+i)] 2 + 4£& im+1 . 



(30) 



In the ansatz of the wavefunction (|6]), the dimensions of the Hilbert space is only 2(N tr + 1). So for each m, we 
only have two eigenvalues for excited states. The other two roots for each m should be omitted. Note that at weak 
coupling, the parity for each eigenstate is fixed and arranged from bottom to above with the order as the first even 
state, then followed by two odd states, two even states, two odd states, and so on. Therefore, the excited states 1 
and 2 are of odd parity, which should be given by Eq 3 ' 4 \ the excited states 3 and 4 are of even parity, then given by 



e[ x ' 2 \ the excited states 5 and 6 are of odd parity, then given by E^ 3 ' 4 \ 
excited states for any m can be summarized as 



and so on. In this way, two eigenvalues for 



E m =m- g 



(-iy 



{D m +l,m+l + D rt 



l )±- x /[i-(-iy 



{D r , 



D 



m+l,m+l ; 



+ 4Ai, m+1 . 



(31) 



Besides, the GS energy is given by Eq 



(i) 



E, 



GS 



1 
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(Di,i + Do,o) 



[1 + (D, 



0,0 



Di. 



(32) 



The FOA results in Eqs. ([3"Tj) and (f3"2"j) have been given directly by the determinant with 2-by-2 block form in Ref. 
(l8| by two present authors and one collaborator previously. We here display the derivation in detail. Especially we 
rule out two pseud solutions for each m by taking the fixed parity of the eigenstates into account. 

Surprisingly Eq. (|3T|) is exactly the same as the previous GRWA result Eq. (20) in Ref. [l6| by Irish. We now 
become aware that the previous GRWA, which were obtained in an alternative way within a lengthy derivation, is 
just FOA in the present scheme. Actually this expression has been derived much earlier within substantially different 
approaches 28] . What is more, we can straightforwardly perform the second-order approximation for the further 
improvement, and extension to the biased case in the present framework, which is however not so easy to operate 
within Irish's approach. To the best of our knowledge, GRWA with finite static bias does not exist in the literature. 

For the finite bias e ^ 0, the parity symmetry is broken with e, with the following notation 



E = x 



9 



-D„ 



(33) 



-D mtTn+ i,w 



-D 



m+l,m-t-l 5 



The determinant can be reduced to 



—x u 
u —x 
v 
v 



v 
■ e v 

1 — x w 
w 1 - 



= 0. 



(34) 
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FIG. 2: The energy levels as a function of coupling constant for for different qubit bias e — 0.1 (upper panel), 0.5(middle 
panel), andl.O (down panel). The values of A are 0.5, 1.0, and 1.5 from left to right column. The present GRWA results (black 
open circles) are compared to the exact (black solid lines), VVP (red dashed lines), and ZOA (blue filled circles) ones. 
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The corresponding quartic equation is 



bx A 



dx + e = 0, 



where 



-2-2e, 
1 + 3e + e 2 



(2v 2 



w 



= (2v 2 



w 2 -e-l)e + 2(u 2 + v 2 ) 



(uw — v 2 ) 2 — u 2 (1 + e) — v 2 e. 



The solutions to this quartic equation are given in the Appendix A. Compared to the exact solutions, we find that 
the second and third roots x 2 and x 3 in Eqs. (IA9|) and (|A10|) are generally the true solutions. The GS energy is given 
by the first root x\ in Eq. (|A8[) for m = 0. We also call the FOA with finite static bias as GRWA. In this way, we can 
calculate the eigenenergies uniquely and straightforwardly, which are shown in Fig. [2] with black circles. It is very 
interesting to find that the present GRWA results are very close to the exact ones in the whole coupling regime for 
wide range of the static bias. Compare to the VVP at static bias e < 1, the present GRWA is obviously much better. 

As stated above, for zero-static bias case e = 0, there is still room to improve by performing the higher order 
approximation. In Ref. fl6j . after a unitary transformation, only the "energy-conserving" one excitation terms like 
their Eq. (18), a generalization of the energy-conserved term in the usual RWA, is kept in their Eqs. (13) and (14), 
so it is called GRWA. Because the present FOA is equivalent to GRWA, so in the second-order approximation, the 
terms beyond their Eq. (18) must be included, so we term this improvement to GRWA as beyond the RWA (BRWA). 
In other words, BRWA can not be implemented within any renormalized RWA form like in Ref. [l6j . 

In the BRWA, the analytical expression can be uniquely and clearly derived within the following procedure. The 
determinant is 



'TO - g 2 - E=f D m , m ) TAn,m + l 
T An : m+1 (TO + 1 - g 2 



•'m,m+l 
m,m+2 



E T D m+ i 



m+1 I 



m+l,m+2 



J2 



[m 



where — (+) for even(odd) parity, which can be simplified as 



-X^u Tx TV 

=Fz (1 - X =F v) 

TV Tz (2-Xtw) 



2 — 9 — E =F D m+ 2. m +2) 



= 0, 



= 0, 



(35) 



where 



E = X + m-g 2 , 

U = D m>m , V = -D m +l,m+l, W = D m+ 2,m+2i 
X = -Dm : m+1, V — D m j yl+ 2; Z = -D m +l,rn+2, 



which gives the following cubic equation 



where 



for even parity, and 



X + bX 2 + cX + d = 0, 



(36) 



[u + v + w) — 3, 

- (x 2 + y 2 + z 2 ) + u (v - 1) + (u + v - 1) (w - 2) . 
(2u — uw + y 2 ) (1 — v) — z 2 u + x 2 (2 — w) + 2xyz, 



— (u + v + w) — 3, 

- (x 2 + y 2 + z 2 ) + u (1 + v) + (it + v + 1) (2 + w) . 
(y 2 — 2u — uw) (1 + v) + z 2 u + x 2 (2 + w) — 2xyz, 
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for odd parity. 

The three different roots to the cubic equation can be found in the Appendix A. In this approximation, we have 
more than one eigenvalues for each eigenstate with fixed parity, which are all true solutions physically, but only some 
of them would be selected. The criterion for the unique formulae for the BRWA is that the solutions are the most close 
to the exact results in the whole parameter regime. In this way, we find that, for even (odd) number m, two roots y\ 
and 1/2 in Eqs. (|A2[) and (|A3|) of the determinant with odd (even) parity would generally give the best eigenvalues for 
the excited states. The GS state is given by the first root y\ of the m = determinant with even parity. Actually, Eq. 
(|35[) has been written out in Ref. [18| by two present authors and one collaborator. But the detailed expression for 
the eigenvalues was not presented. Even the further third approximation was also performed in Ref. [l8[ . The direct 
comparisons between these different order approximations to the GRWA [1(1 HH have not been given, which however 
could reveal the advantage of this scheme. 

We examine the BRWA energy levels against the qubit-oscillator detuning 8 for fixed couplings g — 0.1,0.5, 1.0, 
and 1.5, respectively in Fig. [1] where the GRWA results have been also collected. It is interesting to note that BRWA 
result is always more close to the exact one than GRWA one in all values of the coupling strength, which becomes 
more pronounced with increasing S. 

Due to the counter-rotating wave terms, the eigenfunctions and eigenvalues of the JC model without the RWA 
present an open problem because they are not known in anything like a closed form, even given the exact solutions 
reported recently [22|, H3|. No analytical expressions for the exact eigenvalues are available in the literature, to the 
best of our knowledge. The analytical expressions presented in this paper, which is not exact but work well, might 
be practically useful. 



IV. SUMMARY 



In this paper, by an effective scheme within two displaced bosonic operators with equal positive and negative 
displacements, we study the qubit-oscillator systems analytically in an unified way. Many previous analytical treat- 
ments, such as GRWA, an expansion in the qubit tunneling matrix element in the deep strong coupling regime can 
be recovered in the present scheme. More over, we extend the GRWA to the finite-bias case. The results is much 
better than VVP in the weak and intermediate coupling regime, which is more experimentally interesting. For the 
zero static bias, the GRWA is further improved to BRWA, which is more close to the exact ones at large detuning 
while the GRWA deviates strongly. The analytical expression is explicitly given for future applications. 
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Appendix A: Solutions to univariate cubic and quartic equations 



The univariate cubic equation can be generally expressed as 

x 3 + bx 2 + cx + d = 0. 
Its solutions can be found in any Mathematics manual. If 

r = B 2 - AAC < 0, 

with 

A = b 2 - 3c,B = be- 9d,C = c 2 - 3bd, 

there are three different real roots 



-b-2y/Acos6 , A1 . 

Vi = g » ( A1 ) 

* = (A2) 

„ = -> (« + ¥), (A3) 
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1 (2Ab-W\ 
arccos j= — I . (A4) 



where 

6 = . 

3 V 

The univariate quartic equation can be generally expressed as 

x 4 + bx 3 + cx 2 + dx + e = 0, 
Its four solutions are exactly the four solutions of the following two quadratic equations 

2 i b + z f by — d 



x* + —x+ U + -£_ ) = 0, (A5) 



2 t b — z ( by — d 



x A + ——x +[y- -2 = 0, (A6) 



where z = y/8y + b 2 — 4c and y is the third root 1/3 in Eq. (|A3j) of the following cubic equation 



, c 2 fbd \ eUc-b 2 ) -d 2 
y ~2 y + ) y+ - 8 =°- 



We haver checked that T < in all parameters in the present case. Therefore the four roots are 



1 /, si //■ ,2 I6y(b + z) - 16d 
xi = —^ b + z )--^( b + zf- ^ V -± , (A7) 



1 /, x 1 //, a 16y(6 + z) - 16d /An , 
x 2 = --(b + z) + -^(b + z) 2 - yK -± , (A8) 



1 „ x 1 //, x2 16y (b — z) — 16d , . „. 

^3 = -2(b-z)--^J(b-z) 2 + ^ V ^ , (A9) 



^4 = --(b-z) + -^{b-zf + yy J . (A10) 
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